Dice Check

knitr::opts_chunk$set(echo = TRUE)

pkgs <- c("pwr", "dplyr", "ggplot2", "tidyr", "knitr", "kableExtra", "shiny")
for (i in seq_along(pkgs)) {library(pkgs[i], character.only = TRUE)}

## Warning: package 'pwr' was built under R version 4.4.1

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

## 
## Attaching package: 'kableExtra'

## The following object is masked from 'package:dplyr':
## 
##     group_rows

## Warning: package 'shiny' was built under R version 4.4.1

# DnD character generation option: roll 4d6, take the top 3
# simulating to check distribution of that approach

# function to simulate dice rolls ==============================================
# inputs:
# - counts: vector for number of rolls
# - sides: vector of sides per die (aligned with counts)
# - seed: number to set random number generator seed for replication
#     default (NULL) sets no RNG seed
# - drop.lowest.n: vector to specify if dropping the lowest 'n' results
#     default (NULL) will drop none for any roll
#     a single number will apply that number to any roll
#     a vector of 
# output:
# - list with 1 entry per roll set; each entry has the final rolls
#     (excluding any lowest dropped rolls if specified) and sum of the set
#     the final list entry is the sum across all roll sets
roll <- function(counts, sides, seed = NULL, drop.lowest.n = NULL) {
  if (length(counts) != length(sides)) {stop("counts and sides lengths don't match")}
  if (!all(sides %in% c(4, 6, 8, 10, 100, 12, 20))) {
    stop("sides can only contain 4, 6, 8, 10, 100 (%), 12, 20")}
  set.seed(seed)
  pct_idx <- sides == 100
  
  if (is.null(drop.lowest.n)){d.l.n <- rep(0, length(counts))
  } else {
    if (any(drop.lowest.n >= counts)) {stop("drop.lowest.n must be at most counts - 1")}
    if (length(drop.lowest.n) == 1L) {d.l.n <- rep(drop.lowest.n, length(counts))
    } else {
      if (length(drop.lowest.n) != length(counts)) {
        stop("drop.lowest.n must be one of: same length as counts, a single number, or NULL")
      }
      d.l.n <- drop.lowest.n
      }
  }
  
  # check feasibility of dropping lowest 'n' rolls in each case
  if (any(d.l.n < 0)) {stop("drop.lowest.n cannot be negative")}

  roll_set <- 
    lapply(1:length(sides), function(i) {
      if (pct_idx[i]) {
        ones <- sample(0:9, size = counts[i], replace = TRUE)
        tens <- sample(0:9 * 10, size = counts[i], replace = TRUE)
        rslt <- tens + ones
        rslt[rslt == 0] <- 100
        
        if (d.l.n[i] > 0) {rslt <- (rslt[order(rslt)])[-c(1:d.l.n[i])]}
        rslt <- list("rolls" = rslt, "sum" = sum(rslt))
        rslt
      } else {
        rslt <- sample(1:sides[i], size = counts[i], replace = TRUE)
        if (d.l.n[i] > 0) {rslt <- (rslt[order(rslt)])[-c(1:d.l.n[i])]}
        rslt <- list("rolls" = rslt, "sum" = sum(rslt))
        rslt
        }
      })
  
  names(roll_set) <- paste0(counts, "d", ifelse(pct_idx, "%", sides))
  sums <- 
    vapply(seq_along(roll_set), function(i){roll_set[[i]]$sum}, numeric(1L)) 
  
  roll_set[[length(roll_set) + 1]] <- c("ovr total" = sum(sums))
  
  roll_set
}

# just curious: impact of 3d4 minus 2 versus 1d10 (same min/max)
# data.frame(setup = c("3d4 - 2 vs. 1d10"),
#            sim_seed = 1:1e4) %>%
#   mutate(
#     results = lapply(1:nrow(.), function(i){
#       roll(counts = c(3, 1), sides = c(4, 10), seed = .$sim_seed[i])
#     })
#   ) %>%
#   mutate(
#     res_3d4mns2 = vapply(1:nrow(.), function(i) {
#       .$results[[i]]$`3d4`$sum - 2}, numeric(1L)),
#     res_1d10 = vapply(1:nrow(.), function(i) {
#       .$results[[i]]$`1d10`$sum}, numeric(1L)),
#   ) %>%
#   pivot_longer(cols = starts_with("res_"),
#                names_to = "roll", names_prefix = "res_",
#                values_to = "result") %>%
#   ggplot(aes(x = result, color = roll)) +
#   geom_density() +
#   scale_x_continuous(breaks = 1:10, minor_breaks = FALSE)

# unsurprising - 3d4 - 2 concentrates results closer to center of results range
# (esp. 4-7); 1d10 is more 'swingy'

1 Analysis Plan

• Chi-square Goodness of Fit (GoF) Test

• Bayesian multinomial modeling

2 Dice set notes

• ‘white’: Chessex ‘speckled arctic camo’ 7-die set (CHX 25311); mostly white with grey flecks and black numbering, giving it a ‘granite’ appearance
  • urea
• ‘grey’: Chessex ‘speckled granite’ 7-die set (CHX 25320); white ‘core’ composition with extensive darker grey flecks and some lighter grey flecks, and ‘blood-red’ numbering
  • urea
    
  • retired from use due to the mottled grey/dark red color scheme being hard to read unless under very bright light; re-inking while trying to keep reddish numbering didn’t help
• ‘black’: Chessex ‘speckled urban camo’ 7-die set (CHX 25328); dark grey ‘core’ composition with extensive black flecks and yellow numbering
  • urea
• ‘green’: FanRoll ‘recycled rainbow’ 7-die set; white ‘core’ composition with various-colored flecks throughout and green numbering
  • resin
• ‘teal1’ and ‘teal2’: Gravity Dice individual d20s; teal with silver numbering (same coloring for each; unfortunately before rolling I did not note this, but one has the 15 at the bottom of its facet ‘triangle’ - otherwise they appear identical; based on subsequent test rolls, I believe the ‘15 label at bottom’ die is ‘teal1’)
  • 6061 aluminum
• ‘mushroom’: [Etsy-gifted] individual d20; clear with embedded moss, lavender, and a small mushroom and gold numbering
  • resin
• ‘aubergine’: Die Hard Dice ‘Rerolls’ recycled plastic 7-die set; dark purple (d8 and d20 are closer to ‘plum’, being lighter) with embedded glitter and white numbering
  • acrylic
• ‘ocean’: TheWizardsVault ‘ocean opal’ 7-die set; translucent cyan with gold foil inclusion and gold numbering
  • resin
• ‘ice’: FeverDream Dice ‘ice queen’ 7-die set; opaque white with blue/green microglitter and gold foil inclusion and pastel light blue numbering
  • resin, sharp-edged (all others thus far are round-edged)
    
  • unique among sets thus far, the d4 is a bipyramid (‘shard’) rather than a tetrahedron

3 Rolling Surface

• Initial rolls (white, grey, black, most green rolls): FanRoll velvet dice tray with leather backing (‘rainbow watercolor’); 10” x 10” dimensions, moderate padding

• Later rolls (rest of green, teal[1/2], mushroom, aubergine, ocean, ice rolls): GameGenic Tokens’ Lair magnetic dice tray; 6.5” x 3.75” dimensions, less padding

  • Starting with the ‘ice’ dice, a uniform layer of cork was added to the tray to provide better cushioning for the sharp-edge dice

4 Power Analysis to Determine Sample Size to Vet Dice

sets <- c("white", "grey", "black", "green", "teal1", "teal2", "mushroom", "aubergine", "ocean", "ice")
dice <- paste0("d", c(4, 6, 8, 10, "%", 12, 20))
alpha <- 0.05 # P(reject null[fair die] | actually fair die)
pwr <- 0.80 # P(don't reject null | alt [stated bias exists] is true)

side_levels <- c(0:20, seq(30, 90, by = 10))

# null hypothesis: equal probability for each side
# alternative hypothesis: one side is 'preferred' (100 * alt_pref)% more
#  (by dice structure, that means opposite side dice is 'disliked' by same %)
#  (I'm assuming this holds true for d4 too)
alt_pref <- 0.33333

dice_setup <- 
  data.frame(
    set = rep(sets, each = length(dice)),
    die = rep(dice, times = length(sets))
    ) %>% 
  # only d20 for teal1/teal2/mushroom 'set's
  filter(!((grepl("teal", set) | grepl("mushroom", set)) & 
           die %in% paste0("d", c(4, 6, 8, 10, "%", 12)))) %>% 
  mutate(
    sides_char = ifelse(grepl("\\%", die), 10, sub("d", "", die)),
    sides = as.numeric(sides_char)
  ) %>% 
  select(-sides_char) %>% 
  # note: need separate mutate() calls when the previous mutate-created 
  #       content is needed to create subsequent multi-step entry 
  mutate(
    null = lapply(1:nrow(.), function(i) {
      rep(1 / .$sides[i], times = .$sides[i])
    }),
    alt = lapply(1:nrow(.), function(i) {
      eq_prob <- 1 / .$sides[i]
      prefr <- (1 + alt_pref) * eq_prob
      dislk <- (1 - alt_pref) * eq_prob
      rest <- (1 - (prefr + dislk)) / (.$sides[i] - 2)
      
      c(prefr, dislk, rep(rest, .$sides[i] - 2))
    })
  ) %>% 
  mutate(
    effect_size = vapply(1:nrow(.), function(i) {
      ES.w1(P0 = .$null[[i]], P1 = .$alt[[i]])
    }, numeric(1L))
  ) %>% 
  mutate(
    sample_size = vapply(1:nrow(.), function(i) {
      N_ <- 
      pwr.chisq.test(w = .$effect_size[i], 
                     N = NULL,
                     df = .$sides[i] - 1, 
                     sig.level = alpha,
                     power = pwr)$N
      
     ceiling(N_) # round up
    }, numeric(1L))
  )

5 Checking Dice

5.1 Data Collection

Notes on data collection:

• Rolls must land in the dice tray; any roll outside it will be rejected (this includes rolls hitting the edge of the dice tray and landing outside the tray, rolls landing on the dice tray but then bouncing out, and rolls landing outside the dice tray but then bouncing in; rolls hitting the dice tray side but no other object and landing inside the tray will be included)

• Rolls must be spun in the air; any roll (subjectively) deemed to have been insufficiently ‘spun’ will be rejected (i.e. it felt like I basically just dropped the dice instead of spinning it in the air)

• Dice are rolled one at a time; each result is recorded before picking up the die to conduct the next roll

• If the die doesn’t come to rest fully on one face (i.e. one edge of the die is resting against the side of the dice tray), the die must be (subjectively) deemed completely unambiguous as to which side would be facing up if the tray side were not there for the roll to be recorded

• If the die is rotated while picking it up and before the roll result is noted, the die should be rolled twice without recording the results before resuming recording of roll results (to potentially reduce the risk of bias in the recorded result, should one roll have an association with the prior one)

• The roll technique is described below:
  

• The die starts sitting upright in the open right hand, only slightly cupped with the die generally resting around the base of the middle and ring fingers
  

• The roll is an upward and ‘counterclockwise’ (wrist rotates the pinky to the right) twist, allowing the die to roll from the hand; generally this means the die rolls off the pinky, or occasionally the ring finger
  

• After each roll is recorded, the die is simply picked up and re-rolled; no deliberate care is made to randomize or control the up-facing side of the die (generally it seems to differ from the previous roll’s result due to the picking up action)

# quick plotting function
plot_result <- function(df, fill_color = "dodgerblue2") {
  set_die <- deparse(substitute(df)) # store df name as text
  set_ <- sub("^([a-z]+)_d[0-9|%]{1,2}_df$", "\\1", set_die)
  die_ <- sub("^[a-z]+_(d[0-9|%]{1,2})_df$", "\\1", set_die)
  
  avg_count <- (1 / length(unique(df$result))) * nrow(df)
  
  ggplot(df, aes(x = result)) +
    geom_vline(xintercept = mean(df$result), linetype = "solid") +
    geom_hline(yintercept = avg_count, linetype = "solid") +
    geom_bar(fill = fill_color, color = "black") +
    labs(title = paste("Plot of", set_, die_, "results"),
         subtitle = paste("n = ", nrow(df)),
         caption = paste0("Horizontal line indicates expected count (",
                          sprintf("%.1f", avg_count), 
                          "); Vertical line indicates mean value (",
                          sprintf("%.1f", mean(df$result)),")")) +
    scale_x_continuous(breaks = unique(df$result)) +
    theme_bw()
}


set_fill_map <- 
  c("white" = "ivory3", "grey" = "grey40", 
    "black" = "grey10", "green" = "green4",
    "teal1" = "darkcyan", "teal2" = "steelblue4",
    "mushroom" = "gold3", "aubergine" = "#614051",
    "ocean" = "cyan2", "ice" = "powderblue")

5.1.1 White Dice

5.1.2 Grey Dice

5.1.3 Black Dice

5.1.4 Green Dice

5.1.5 Individual d20s

5.1.6 Aubergine Dice

5.1.7 Ocean Dice

5.1.8 Ice Dice

6 Aggregated Analyses

# loading the above dice-roll data
dice_roll_df <- 
  read.csv(file = "dice_rolls.csv")

# function to identify streaks in the data
find_run_lengths <- function(x, streaks.only = TRUE){
  run_len_encoding <- rle(x)
  rle_df <- data.frame(value = run_len_encoding$values,
                       length = run_len_encoding$lengths)

  if (streaks.only) {
    rle_df[rle_df$length > 1L,]
  } else {rle_df}
}

# function for multinomial proportion modelling
rdirichlet <- # from the LearnBayes package
  function (n, par) {
    k = length(par)
    z = array(0, dim = c(n, k))
    s = array(0, dim = c(n, 1))
    
    for (i in 1:k) {
      z[, i] = rgamma(n, shape = par[i])
      s = s + z[, i]
    }
    
    for (i in 1:k) {
      z[, i] = z[, i] / s
    }
    
    return(z)
  }

dice_roll_df_ <- 
  dice_roll_df %>% 
  mutate(set = factor(set, levels = sets),
         die = factor(die, levels = paste0("d", c(4, 6, 8, 10, "%", 12, 20)))
  ) %>% 
  nest_by(set, die) %>% 
  mutate(
    die_val = sub("d", "", die),
    die_val = ifelse(die_val == "%", 100, as.numeric(die_val)),
    n_obs = length(data$result),
    expected_count_per_side = 
      case_when(die_val == 100 ~ n_obs / 10,
                TRUE ~ n_obs / die_val),
    mean_ = mean(data$result),
    expected_mean = 
      case_when(die_val == 100 ~ mean(0:9 * 10),
                die_val == 10 ~ mean(1:10),
                TRUE ~ mean(1:die_val)),
    sd_ = sd(data$result),
    expected_sd = 
      case_when(die_val == 100 ~ sd(0:9 * 10),
                die_val == 10 ~ sd(1:10),
                TRUE ~ sd(1:die_val)),
    result_tbl = list(table(data$result)),
    chisq_pval = chisq.test(result_tbl)$p.value,
    chisq_stat_sig_0.05 = chisq_pval <= 0.05,
    roll_streaks = list(find_run_lengths(data$result, streaks.only = TRUE)),
    roll_streak_tbl = list(table(roll_streaks))
    )

6.1 ‘Fairest’ and ‘Least Fair’

# checking which dice per side across sets have largest p-values
fairest_dice_top4 <- 
  dice_roll_df_ %>% 
  group_by(die) %>% 
  arrange(desc(chisq_pval)) %>% 
  mutate(p_val_rank = 1:n()) %>% 
  ungroup() %>% 
  filter(p_val_rank %in% c(1:4)) %>% 
  select(set, die, chisq_pval, p_val_rank) %>% 
  arrange(die, p_val_rank, set)

fairest_dice_bot3 <- 
  dice_roll_df_ %>% 
  group_by(die) %>% 
  arrange(desc(chisq_pval)) %>% 
  mutate(p_val_rank = 1:n(), max_rank = n()) %>% 
  ungroup() %>% 
  mutate(bad = p_val_rank == max_rank - 2,
         worse = p_val_rank == max_rank - 1,
         worst = p_val_rank == max_rank) |> 
  filter(bad | worse | worst) %>% 
  select(set, die, chisq_pval, p_val_rank) %>% 
  arrange(die, desc(p_val_rank), set)

# table formatting
format_table <- function(x){
  x |> 
  kbl(format = "html") |>
  kable_styling(full_width = FALSE, 
                bootstrap_options = c("striped", "bordered"))
}

noquote("Fairest Dice - Top 4")

## [1] Fairest Dice - Top 4

fairest_dice_top4 |> 
  mutate(chisq_pval = round(chisq_pval, 3)) |> 
  format_table()

set	die	chisq_pval	p_val_rank
grey	d4	0.900	1
aubergine	d4	0.881	2
ocean	d4	0.641	3
black	d4	0.632	4
aubergine	d6	0.802	1
green	d6	0.612	2
grey	d6	0.561	3
black	d6	0.253	4
aubergine	d8	0.982	1
grey	d8	0.910	2
white	d8	0.411	3
black	d8	0.386	4
grey	d10	0.829	1
white	d10	0.804	2
black	d10	0.793	3
green	d10	0.556	4
grey	d%	0.624	1
white	d%	0.363	2
black	d%	0.305	3
aubergine	d%	0.130	4
white	d12	0.507	1
green	d12	0.393	2
ice	d12	0.387	3
black	d12	0.321	4
aubergine	d20	0.097	1
ice	d20	0.002	2
grey	d20	0.001	3
teal2	d20	0.000	4

noquote("Fairest Dice - Sets by # Appearances and Average Ranking in Top 4")

## [1] Fairest Dice - Sets by # Appearances and Average Ranking in Top 4

fairest_dice_top4 |> 
  group_by(set) |> 
  summarize(top4_ct = n(), 
            avg_rank = mean(p_val_rank) |> round(3), 
            .groups = "drop") |> 
  arrange(desc(top4_ct), avg_rank) |> 
  format_table()

set	top4_ct	avg_rank
grey	6	1.833
black	6	3.667
aubergine	5	1.800
white	4	2.000
green	3	2.667
ice	2	2.500
ocean	1	3.000
teal2	1	4.000

# grey, then aubergine, then black, then white, then green 
# (descending order of 'fairness')
# aside from d20s, no real concerns of fairness for any observations
# aubergine d20 is astonishingly 'fair'!
# surprised the ice d20 is second fairest...

noquote("Least Fair Dice - Bottom 3")

## [1] Least Fair Dice - Bottom 3

fairest_dice_bot3 |> 
  mutate(chisq_pval = round(chisq_pval, 3)) |>
  format_table()

set	die	chisq_pval	p_val_rank
ice	d4	0.588	7
white	d4	0.597	6
green	d4	0.614	5
ocean	d6	0.024	7
white	d6	0.068	6
ice	d6	0.187	5
green	d8	0.040	7
ice	d8	0.170	6
ocean	d8	0.351	5
ocean	d10	0.051	7
ice	d10	0.070	6
aubergine	d10	0.133	5
ice	d%	0.004	7
green	d%	0.055	6
ocean	d%	0.105	5
ocean	d12	0.014	7
aubergine	d12	0.079	6
grey	d12	0.287	5
green	d20	0.000	10
white	d20	0.000	9
black	d20	0.000	8

noquote("Least Fair Dice, Arranged by p-value")

## [1] Least Fair Dice, Arranged by p-value

fairest_dice_bot3 |> 
  arrange(chisq_pval) |> 
  mutate(chisq_pval = round(chisq_pval, 3)) |>
  format_table()

set	die	chisq_pval	p_val_rank
green	d20	0.000	10
white	d20	0.000	9
black	d20	0.000	8
ice	d%	0.004	7
ocean	d12	0.014	7
ocean	d6	0.024	7
green	d8	0.040	7
ocean	d10	0.051	7
green	d%	0.055	6
white	d6	0.068	6
ice	d10	0.070	6
aubergine	d12	0.079	6
ocean	d%	0.105	5
aubergine	d10	0.133	5
ice	d8	0.170	6
ice	d6	0.187	5
grey	d12	0.287	5
ocean	d8	0.351	5
ice	d4	0.588	7
white	d4	0.597	6
green	d4	0.614	5

# green d20 is clearly not random, nor is white d20
# for other-sided dice even the 'worst' performers are solid except 
# 'ice' d10 slightly a bit biased against '0', EV is ok (5.4 vs. 'fair' 5.5)
# 'ice' d% is very unlikely to roll 00, otherwise looks great
# 'ocean' d12 ... expected value (6.3) isn't wildly off 'random' EV (6.5)
# 'ocean' d6  ... expected value (3.5) exactly matches 'random' EV, 
#   just more 3s than 1s 

6.2 Mean and SD of Observed Rolls

# observed vs. expected mean / SD
dice_mean_df_tbl <- 
  dice_roll_df_ %>% 
  select(set, die, observed_mean = mean_, expected_mean) %>% 
  mutate("pct_diff" = (observed_mean - expected_mean) / expected_mean * 100) %>% 
  pivot_longer(cols = ends_with("_mean"), names_to = "mean", values_to = "value") %>% 
  mutate(type = sub("_mean", "", mean), var = "mean") %>% 
  select(-mean)

dice_sd_df_tbl <- 
  dice_roll_df_ %>% 
  select(set, die, observed_sd = sd_, expected_sd) %>% 
  mutate("pct_diff" = (observed_sd - expected_sd) / expected_sd * 100) %>% 
  pivot_longer(cols = ends_with("_sd"), names_to = "sd", values_to = "value") %>% 
  mutate(type = sub("_sd", "", sd), var = "sd") %>% 
  select(-sd)

dice_mean_sd_tbl <- 
  full_join(dice_mean_df_tbl, dice_sd_df_tbl, 
            by = c("set", "die", "var", "type", "value", "pct_diff")) %>% 
  mutate(gt_10pct_diff = abs(pct_diff) > 10)

mean_sd_by_set_plots <- 
  lapply(seq_along(sets), function(i) {
    ggplot(data = dice_mean_sd_tbl[dice_mean_sd_tbl$set == sets[i],],
           aes(x = type, ymin = 0, ymax = value, y = value, color = gt_10pct_diff)) +
      facet_grid(rows = vars(die), cols = vars(var), scales = "free_y") +
      geom_pointrange() +
      labs(title = "Expected vs. Observed Mean/SD of Dice Rolls",
           subtitle = paste(sets[i], "dice"),
           color = ">10% diff.\nvs. expected") +
      scale_color_manual(values = c("FALSE" = "grey30", "TRUE" = "orange3")) +
      theme_bw()
  })

# focus on cases of interest
focus_mean_sd_plots <- 
  dice_mean_sd_tbl %>% 
  filter(gt_10pct_diff) %>% 
  nest_by(set, die) %>% 
  mutate(
    plot_ = list(
      ggplot(data = data,
           aes(x = type, ymin = 0, ymax = value, y = value)) +
      facet_wrap(facets = vars(var), scales = "free_y", ncol = 2) +
      geom_pointrange(color = "orange3") +
      labs(title = "Expected vs. Observed Mean/SD of Dice Rolls",
           subtitle = paste(set, die)) +
        theme_bw()
    )
  ) %>% 
  pull(plot_)

for (i in seq_along(mean_sd_by_set_plots)) {
  if (i == 1) {cat("All Mean/SD Plots\n\n")}
  cat("\n\n")
  print(mean_sd_by_set_plots[[i]])
}

All Mean/SD Plots

for (i in seq_along(focus_mean_sd_plots)) {
  if (i == 1) {cat("10+% difference Mean/SD Plots\n\n")}
  cat("\n\n")
  print(focus_mean_sd_plots[[i]])
}

10+% difference Mean/SD Plots

6.3 Roll Streaks

streak_color_map <- 
  c("1" = "blue4", "2" = "blue2", "3" = "dodgerblue3",
    "4" = "firebrick1", "5" = "firebrick", "6" = "red4")

dice_streaks_plots <- 
  lapply(1:nrow(dice_roll_df_), function(i) {
    df_0 <- dice_roll_df_$roll_streak_tbl[[i]]
    df_ <- data.frame(value = rep(rownames(df_0), times = ncol(df_0)),
                      strk_len = rep(colnames(df_0), each = nrow(df_0)))
    df_$value <- as.factor(as.integer(df_$value))
    
    df_$count <- as.integer(df_0)

    ggplot(df_, aes(x = value, y = count, fill = strk_len)) +
      geom_col(position = position_dodge(width = 0.7)) +
      labs(title = "Streak length by value",
           subtitle = paste(dice_roll_df_$set[i], dice_roll_df_$die[i])) +
      scale_fill_manual(values = streak_color_map) +
      theme_bw()
  })

streak_sim <- 
  lapply(1:100, function(i) {
    roll_sim <- roll(counts = 1851, sides = 20, seed = i)[[1]]$rolls
    rle_ <- find_run_lengths(roll_sim, streaks.only = TRUE)
    rle_
  }) %>% 
  bind_rows() %>% 
  table()

streak_sim_df <- 
  data.frame(value = rep(rownames(streak_sim), times = ncol(streak_sim)),
             strk_len = rep(colnames(streak_sim), each = nrow(streak_sim)))
streak_sim_df$value <- as.factor(as.integer(streak_sim_df$value))
streak_sim_df$count <- as.integer(streak_sim)

streak_sim_plot <- 
  ggplot(streak_sim_df, aes(x = value, y = count, fill = strk_len)) +
  geom_col(position = position_dodge(width = 0.7)) +
  labs(title = "Streak length by value",
       subtitle = "Simulated d20 rolls",
       caption = "1,851 rolls simulated 100 times") +
  scale_fill_manual(values = streak_color_map) +
  theme_bw()

dice_streak_focus <- c(7, 14, 21, 28, 29, 30, 31, 38, 40, 45, 52)

focused_streak_plots <- 
  lapply(seq_along(dice_streak_focus), function(i){
    dice_streaks_plots[[dice_streak_focus[i]]]
    })

cat("\n\nSimulated d20 rolls for comparison\n\n")

Simulated d20 rolls for comparison

print(streak_sim_plot)

cat("\n\nSelected Streak Plots\n\n")

Selected Streak Plots

for (i in seq_along(focused_streak_plots)) {
  cat("\n\n")
  print(focused_streak_plots[[i]])
}

6.4 Observed vs. Expected for Stat. Sig. Dice

Note: The ‘aubergine’ d20 is (amazingly) not statistically significant; it’s included for comparison against the other d20s, all of which were stat. sig. in the chi-squared test.

# plots of (Chi square HoV) stat-sig dice - side-by-side for same-die sets
comparison_plots_df <- 
  dice_roll_df_ %>%
  filter(chisq_stat_sig_0.05 |
         (set == "aubergine" & die == "d20")) %>% 
  nest_by(die)

set_color_map <- 
  c("white" = "black", "grey" = NA_character_,  
    "black" = NA_character_, "green" = NA_character_, 
    "teal1" = NA_character_, "teal2" = NA_character_,
    "mushroom" = NA_character_, "aubergine" = NA_character_,
    "ocean" = NA_character_, "ice" = NA_character_)

comparison_plots <- 
  lapply(1:nrow(comparison_plots_df), function(i) {
    df_ <- comparison_plots_df[i,]
    summary_data <- df_$data
    results_data <- 
      lapply(summary_data, `[[`, "data") %>% 
      bind_rows() %>% 
      mutate(set = rep(summary_data[[1]]$set, times = summary_data[[1]]$n_obs))
    
    x_vals <- unique(results_data$result)
    x_vals <- x_vals[order(x_vals)]
    
    ggplot(results_data, aes(x = result, fill = set, color = set)) +
      geom_hline(data = summary_data[[1]], 
                 aes(yintercept = expected_count_per_side),
                 linetype = "dashed") +
      geom_bar(position = position_dodge(width = 0.7)) +
      labs(title = paste("Observed Counts for", df_$die, "Rolls by Set"),
           subtitle = "Horizontal line: expected roll count (fair die)") +
      scale_x_continuous(breaks = x_vals, minor_breaks = FALSE) +
      scale_color_manual(values = set_color_map, guide = "none") +
      scale_fill_manual(values = set_fill_map) +
      theme_bw()
  })

for (i in 1:(length(comparison_plots) - 1)) {
  cat("\n\n")
  print(comparison_plots[[i]])
}

cat("\n\n")

print(comparison_plots[[length(comparison_plots)]])

6.5 Checking Stat. Sig. d20 Results If Collecting Smaller Samples

# simulating p-values from subsampling observed rolls - checking if
#  stat. sig. p-vals are more a result of large sample size
sample_chisq_pval <- 
  function(data, subset_var, value_var, nsim = 1e4, n = 500, seed = 42) {
    p_val_vec <- vector("numeric", nsim)
    
    data_sets <- lapply(1:length(unique(data[[subset_var]])), function(i) {
      set <- data[data[[subset_var]] == unique(data[[subset_var]])[i], "set"][1]
      data <- data[data[[subset_var]] == unique(data[[subset_var]])[i], 
                   value_var]
      list("set" = set, "data" = data)
    })
    
    p_val_mat <- 
      matrix(nrow = nsim, ncol = length(data_sets) + 1,
             dimnames = list(NULL, 
                             c("seed", 
                               vapply(data_sets, `[[`, character(1L), 1)))
             )
    
    for (i in 1:nrow(p_val_mat)) {
      for (j in 2:ncol(p_val_mat)) {
        set.seed(seed + i)
        p_val_mat[i, 1] <- seed + i 
        p_val_mat[i, j] <- 
          chisq.test(
            table(
              sample(data_sets[[j - 1]]$data, size = n, replace = FALSE)
              )
            )$p.value
      }
    }
    p_val_mat
  }


d20_sample_sim_n1000 <- 
  sample_chisq_pval(
    data = dice_roll_df[dice_roll_df$die == "d20", c("set", "result")],
    subset_var = "set",
    value_var = "result",
    nsim = 5e3, n = 1000, seed = 42
  )

d20_sample_sim_n750 <- 
  sample_chisq_pval(
    data = dice_roll_df[dice_roll_df$die == "d20", c("set", "result")],
    subset_var = "set",
    value_var = "result",
    nsim = 5e3, n = 750, seed = 42
  )

d20_sample_sim_n500 <- 
  sample_chisq_pval(
    data = dice_roll_df[dice_roll_df$die == "d20", c("set", "result")],
    subset_var = "set",
    value_var = "result",
    nsim = 5e3, n = 500, seed = 42
  )

d20_sample_sim_n1000_df <- 
  apply(d20_sample_sim_n1000, 2, quantile, probs = c(0.10, 0.50, 0.90)) %>% 
  as.data.frame() %>% 
  select(-seed) %>% 
  mutate(quantile = rownames(.)) %>% 
  pivot_longer(cols = all_of(sets), names_to = "set", values_to = "value") %>% 
  mutate(n = 1000) %>% 
  pivot_wider(names_from = "quantile", values_from = "value")

d20_sample_sim_n750_df <- 
  apply(d20_sample_sim_n750, 2, quantile, probs = c(0.10, 0.50, 0.90)) %>% 
  as.data.frame() %>% 
  select(-seed) %>% 
  mutate(quantile = rownames(.)) %>% 
  pivot_longer(cols = all_of(sets), names_to = "set", values_to = "value") %>% 
  mutate(n = 750) %>% 
  pivot_wider(names_from = "quantile", values_from = "value")

d20_sample_sim_n500_df <- 
  apply(d20_sample_sim_n500, 2, quantile, probs = c(0.10, 0.50, 0.90)) %>% 
  as.data.frame() %>% 
  select(-seed) %>% 
  mutate(quantile = rownames(.)) %>% 
  pivot_longer(cols = all_of(sets), names_to = "set", values_to = "value") %>% 
  mutate(n = 500) %>% 
  pivot_wider(names_from = "quantile", values_from = "value")

rm(d20_sample_sim_n1000, d20_sample_sim_n750, d20_sample_sim_n500)

merge_df <- function(df_list) {
    Reduce(function(x, y) merge(x, y, all = TRUE), df_list)
    }

d20_sample_sim_df <- 
  merge_df(list(d20_sample_sim_n500_df, 
                d20_sample_sim_n750_df,
                d20_sample_sim_n1000_df)) 

write.csv(x = d20_sample_sim_df,
          file = "d20_sample_sim.csv", 
          row.names = FALSE)

d20_sample_sim_df <- 
  read.csv("d20_sample_sim.csv")

set_d20_ranks <- 
  d20_sample_sim_df |> 
  group_by(set) |> 
  summarize(X50. = mean(X50.), .groups = "drop") |> 
  arrange(desc(X50.)) |> 
  pull(set)

set_d20_ranks_abbrevmap <- 
  c("aubergine" = "aub",
    "ice" = "ice",
    "grey" = "gry",
    "teal2" = "tl2",
    "ocean" = "ocn",
    "mushroom" = "msh",
    "teal1" = "tl1",
    "black" = "blk",
    "white" = "wht",
    "green" = "grn")

d20_sample_sim_df |> 
  mutate(set = factor(set, levels = set_d20_ranks)) |> 
  ggplot(aes(x = set, ymin = X10., y = X50., ymax = X90., color = set)) +
  geom_hline(yintercept = 0.05, linetype = "dashed") +
  geom_pointrange() +
  facet_wrap(facets = vars(n), ncol = 2) +
  labs(title = "p-values of Chi-square test of d20 by set",
       subtitle = "Summary of 1,000 samples of size 500 / 1000",
       y = "Chi Square p-value",
       caption = "Sampling without replacement; 10/50/90% p-value summary; ref. line at 0.05") +
  scale_color_manual(values = set_fill_map, guide = "none") +
  theme_bw() + 
  theme(panel.grid.major.x = element_blank(),
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1))

6.6 Bayesian Multinomial Analysis (d20)

6.6.1 Analysis and Posterior Probability Density

# brms note for Dirichlet model:
#https://discourse.mc-stan.org/t/dirichlet-regresion-using-brms/8591

set.seed(42)
d20_bayesian_pcts_df <- 
  data.frame(
    set = sets,
    die = "d20"
  ) %>% 
  left_join(
    y = dice_roll_df_ %>% 
      filter(die == "d20") %>% 
      select(set, die, results = result_tbl),
    by = c("set", "die")
    ) %>% 
  mutate(
    thetas_df = lapply(1:nrow(.), function(i) {
      set_ <- .$set[i]
      die_ <- .$die[i]
      
      sim_df <- 
        rdirichlet(n = 1e3, par = unlist(unname(.$results[i])) + 1) %>% 
        as.data.frame()
      
      colnames(sim_df) <- paste0("_", 1:20)
      sim_df %>% 
        pivot_longer(cols = everything(), names_to = "side", values_to = "value",
                     names_prefix = "_") %>% 
        mutate(set = set_, die = die_)
    })
  )

d20_bayesian_pcts_df <- 
  lapply(1:nrow(d20_bayesian_pcts_df), function(i) {
    d20_bayesian_pcts_df$thetas_df[[i]]
    }) %>% 
  bind_rows() %>% 
  mutate(side = factor(side, levels = as.character(1:20)))

d20_bayesian_pcts_df$iter <- 
  rep(
    rep(1:(nrow(d20_bayesian_pcts_df) / 
             (20 * length(unique(d20_bayesian_pcts_df$set)))), 
      each = 20),
    times = length(unique(d20_bayesian_pcts_df$set))
  )
 
set_lntyp_map <- 
  c("white" = "solid", "grey" = "solid", 
    "black" = "solid", "green" = "solid",
    "teal1" = "solid", "teal2" = "solid",
    "mushroom" = "solid", "aubergine" = "dashed",
    "ocean" = "solid", "ice" = "solid")
   
d20_bayesian_pcts_plot_top5 <- 
  d20_bayesian_pcts_df |> 
  filter(set %in% set_d20_ranks[1:5]) |> 
  ggplot(aes(x = value, color = set, linetype = set)) +
  facet_wrap(facets = vars(side), ncol = 4) +
  geom_vline(xintercept = 1 / 20, linetype = "dashed") +
  geom_density(linewidth = 0.8) +
  labs(title = "Bayesian posterior density of d20 roll probabilities by die set (top 5)",
       subtitle = "uniform Dirichlet prior, 1851 rolls per die",
       caption = "Reference line at 5%") +
  scale_x_continuous(labels = scales::percent) +
  scale_color_manual(values = set_fill_map) +
  scale_linetype_manual(values = set_lntyp_map) +
  theme_bw()

d20_bayesian_pcts_plot_nottop5 <- 
  d20_bayesian_pcts_df |> 
  filter(!(set %in% set_d20_ranks[1:5])) |> 
  ggplot(aes(x = value, color = set, linetype = set)) +
  facet_wrap(facets = vars(side), ncol = 4) +
  geom_vline(xintercept = 1 / 20, linetype = "dashed") +
  geom_density(linewidth = 0.8) +
  labs(title = "Bayesian posterior density of d20 roll probabilities by die set (not top 5)",
       subtitle = "uniform Dirichlet prior, 1851 rolls per die",
       caption = "Reference line at 5%") +
  scale_x_continuous(labels = scales::percent) +
  scale_color_manual(values = set_fill_map) +
  scale_linetype_manual(values = set_lntyp_map) +
  theme_bw()

cat("\nTop 5 d20 sets\n")

Top 5 d20 sets

print(d20_bayesian_pcts_plot_top5)

cat("\n\nRemaining d20 sets\n")

Remaining d20 sets

print(d20_bayesian_pcts_plot_nottop5)

6.6.2 Credible Intervals from Posterior Probabilities

d20_bayesian_credint_df <- 
  d20_bayesian_pcts_df |> 
  group_by(set, side) |> 
  summarize(
    credint_pct2.5 = quantile(value, 0.025),
    credint_pct50 = median(value),
    credint_pct97.5 = quantile(value, 0.975),
    .groups = "drop"
  )

d20_bayesian_credint_df <- 
  mutate(d20_bayesian_credint_df, 
         set_abbrev = vapply(1:nrow(d20_bayesian_credint_df), function(i) {
           set_d20_ranks_abbrevmap[grep(d20_bayesian_credint_df$set[i],
                                        names(set_d20_ranks_abbrevmap))]
           }, character(1L))
    ) |> 
  mutate(set_abbrev = factor(set_abbrev, levels = unname(set_d20_ranks_abbrevmap)),
         set = factor(set, levels = names(set_d20_ranks_abbrevmap)))

# plotting the credible intervals
d20_bayesian_credint_df |> 
  ggplot(aes(x = set_abbrev, 
             ymin = credint_pct2.5, y = credint_pct50, ymax = credint_pct97.5,
             color = set)) + 
  facet_wrap(facets = vars(side), ncol = 4) +
  geom_hline(yintercept = 1 / 20, linetype = "dashed") +
  geom_pointrange(linewidth = 1) +
  labs(title = "d20 95% Credible Intervals with Median by Set and Side",
       y = "Posterior Probability (%)",
       caption = "Reference line at 5%") +
  scale_y_continuous(labels = scales::percent) +
  scale_color_manual(values = set_fill_map) +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 40, hjust = 1, vjust = 1))

6.6.3 Credible Interval Zoom-in: Rolling 1 vs 20, and ‘Low’ vs. ‘High’ Results

# deep dive: posterior probability 1 vs. 20
d20_bayesian_credint_df |> 
  filter(side %in% c("1", "20")) |>
  ggplot(aes(x = set, 
             ymin = credint_pct2.5, y = credint_pct50, ymax = credint_pct97.5,
             color = set, group = side)) +
  geom_hline(yintercept = 1 / 20, linetype = "dashed") +
  geom_pointrange(fatten = 10, shape = 22, linewidth = 1, fill = "white", 
                  position = position_dodge(width = 0.5)) +
  geom_text(aes(label = side), size = 3, color = "black",
            position = position_dodge(width = 0.5), vjust = 0.5) +
  labs(title = "d20 95% Credible Intervals with Median by Set and Side",
       subtitle = "Zoom in for P(1) and P(20)",
       y = "Posterior Probability (%)",
       caption = "Reference line at 5%") +
  scale_y_continuous(labels = scales::percent) +
  scale_color_manual(values = set_fill_map, guide = "none") +
  theme_bw()

# probability 1-3 vs 18-20
btm_top <- 3

d20_bayesian_btm_df <- 
  d20_bayesian_pcts_df |> 
  filter(side %in% c(1:btm_top)) |> 
  group_by(set, iter) |> 
  summarize(value = sum(value), .groups = "drop") |> 
  mutate(grp = paste("1 to", btm_top)) |> 
  group_by(grp, set) |> 
  summarize(
    credint_pct2.5 = quantile(value, 0.025),
    credint_pct50 = median(value),
    credint_pct97.5 = quantile(value, 0.975),
    .groups = "drop"
  )

d20_bayesian_top_df <- 
  d20_bayesian_pcts_df |> 
  filter(side %in% c((20 - btm_top + 1):20)) |> 
  group_by(set, die, iter) |> 
  summarize(value = sum(value), .groups = "drop") |> 
  mutate(grp = paste((20 - btm_top + 1), "to 20")) |> 
  group_by(grp, set) |> 
  summarize(
    credint_pct2.5 = quantile(value, 0.025),
    credint_pct50 = median(value),
    credint_pct97.5 = quantile(value, 0.975),
    .groups = "drop"
  )

d20_bayesian_lohi_df <-
  full_join(d20_bayesian_btm_df, d20_bayesian_top_df,
            by = intersect(colnames(d20_bayesian_btm_df),
                           colnames(d20_bayesian_top_df)))
d20_bayesian_lohi_df |> 
  mutate(set = factor(set, levels = names(set_d20_ranks_abbrevmap))) |> 
  ggplot(aes(x = set, 
             ymin = credint_pct2.5, y = credint_pct50, ymax = credint_pct97.5,
             color = set, shape = grp)) +
  geom_hline(yintercept = btm_top / 20, linetype = "dashed") +
  geom_pointrange(fatten = 5, linewidth = 1, fill = "white", 
                  position = position_dodge(width = 0.5)) +
  labs(title = "d20 95% Credible Intervals with Median by Set and Side",
       subtitle = paste("Zoom in for lowest and highest", btm_top, "sides"),
       y = "Posterior Probability (%)",
       caption = paste("Reference line at", paste0(btm_top / 20 * 100, "%")),
       shape = NULL) +
  scale_y_continuous(labels = scales::percent) +
  scale_color_manual(values = set_fill_map, guide = "none") +
  scale_shape_manual(values = c(25, 24)) +
  theme_bw()

6.6.4 Median Posterior Probability of (d20 side) by Set

Dynamic Shiny app

Note: this was an experiment in embedding a Shiny app in an R Markdown file; the HTML file needs to be rendered locally from the .Rmd file for this to work. IF DOING THIS, ADD ‘runtime: shiny’ AT THE LAST LINE OF THE .RMD YAML HEADER

# Shiny app radar plot for median posterior probability of sides per d20
shinyApp(

  ui = fluidPage(
    selectInput("set", "Dice set(s) to emphasize:",
                choices = unique(d20_bayesian_pcts_df$set)[order(unique(d20_bayesian_pcts_df$set))],
                selected = unique(d20_bayesian_pcts_df$set)[order(unique(d20_bayesian_pcts_df$set))][1],
                multiple = TRUE),
    checkboxInput("arrange_desc", "Arrange sides by median posterior probability"),
    plotOutput("radar_plot", width = "100%")
  ),

  server = function(input, output) {
    median_prob_df <- 
      d20_bayesian_pcts_df |> 
      group_by(die, set, side) |> 
      summarize(median_post_prob = median(value), .groups = "drop")
    
    sorted_side_order <- 
      reactive({
        if (input$arrange_desc) {
          median_prob_df |> 
            filter(set %in% input$set) |>
            group_by(side) |> 
            summarize(mean_medn_prob = mean(median_post_prob), 
                      .groups = "drop") |> 
            mutate(side = as.integer(side)) |> 
            arrange(desc(mean_medn_prob)) |> 
            pull(side)
          } else { 1:20 }
      })
    
    medn_prob_df <- 
      reactive({
        mutate(median_prob_df, side = factor(side, levels = sorted_side_order()))
      })
      
    emph_prob_df <-  
      reactive({ 
        median_prob_df |> 
        filter(set %in% input$set) |> 
        mutate(side = factor(side, levels = sorted_side_order()))
        }) 
    
    output$radar_plot = renderPlot({
        ggplot(data = emph_prob_df(), aes(x = side, y = median_post_prob, color = set, group = set)) +
        geom_hline(yintercept = 1 / 20, linetype = "dashed") +
        geom_line(data = medn_prob_df(), color = "grey") +
        geom_point(data = medn_prob_df(), color = "grey") +
        geom_line(data = emph_prob_df(), linewidth = 1) +
        geom_point(data = emph_prob_df(), size = 3) +
        labs(title = "Radar-like Plot of Median Posterior Probability of (side)",
             subtitle = "All d20 in grey, selected set(s) in color",
             x = NULL, y = NULL, color = NULL,
             caption = "Reference line at 5%") +
        scale_y_continuous(labels = scales::percent) +
        scale_color_manual(values = set_fill_map) +
        theme_bw() +
        theme(axis.text.x = element_text(size = 10),
              axis.text.y = element_text(size = 10),
              plot.caption = element_text(size = 10),
              legend.text = element_text(size = 10)) +
        coord_polar()
    })
  },

  options = list(height = 550)
)

Shiny applications not supported in static R Markdown documents

Static Plots

mdn_post_prob_plts <- 
  lapply(seq_along(set_d20_ranks), function(i) {
    median_prob_df <- 
      d20_bayesian_pcts_df |> 
      group_by(die, set, side) |> 
      summarize(median_post_prob = median(value), .groups = "drop")
    
    sorted_side_order <- 
      median_prob_df |> 
      filter(set == set_d20_ranks[i]) |>
      group_by(side) |> 
      summarize(mean_medn_prob = mean(median_post_prob), 
                .groups = "drop") |> 
      mutate(side = as.integer(side)) |> 
      arrange(desc(mean_medn_prob)) |> 
      pull(side)
    
    medn_prob_df <- 
      mutate(median_prob_df, side = factor(side, levels = sorted_side_order))
      
    emph_prob_df <-  
      median_prob_df |> 
      filter(set == set_d20_ranks[i]) |> 
      mutate(side = factor(side, levels = sorted_side_order))
    
    ggplot(data = emph_prob_df, aes(x = side, y = median_post_prob, color = set, group = set)) +
        geom_hline(yintercept = 1 / 20, linetype = "dashed") +
        geom_line(data = medn_prob_df, color = "grey") +
        geom_point(data = medn_prob_df, color = "grey") +
        geom_line(data = emph_prob_df, linewidth = 1) +
        geom_point(data = emph_prob_df, size = 3) +
        labs(title = "Radar-like Plot of Median Posterior Probability of (side)",
             subtitle = "All d20 in grey, selected set in color",
             x = NULL, y = NULL, color = NULL,
             caption = "Reference line at 5%") +
        scale_y_continuous(labels = scales::percent) +
        scale_color_manual(values = set_fill_map) +
        theme_bw() +
        coord_polar()
  })

for (i in seq_along(mdn_post_prob_plts)) {
  print(mdn_post_prob_plts[[i]])
}